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Abstract 

A bivariate version of the multicanonical Monte Carlo method and its apphcation to the simula- 
tion of the three-dimensional it J Ising spin glass are described. We found the autocorrelation time 
associated with this particular multicanonical method was approximately proportional to the sys- 
tem volume, which is a great improvement over previous methods applied to spin-glass simulations. 
The principal advantage of this version of the multicanonical method, however, was its ability to 
access information predictive of low-temperature behavior. At low temperatures we found results 
on the three-dimensional it J Ising spin glass consistent with a double degeneracy of the ground- 
state: the order-parameter distribution function P{q) converged to two delta-function peaks and 
the Binder parameter approached unity as the system size was increased. With the same density 
of states used to compute these properties at low temperature, we found their behavior changing 
as the temperature is increased towards the spin glass transition temperature. Just below this 
temperature, the behavior is consistent with the standard mean-field picture that has an infinitely 
degenerate ground state. Using the concept of zero-energy droplets, we also discuss the structure 
of the ground-state degeneracy. The size distribution of the zero-energy droplets was found to 
produce the two delta-function peaks of P{q)- 
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I. INTRODUCTION 



Although a quarter century has passed since the pioneering work on spin glasses by Ed- 
wards and Anderson , many fundamental problems remain unsolved even for the simplest 
models of these materials 0. In the present paper, we report the results of Monte Carlo 
simulations of the three-dimensional ±J Ising spin glass, focusing on the nature of the 
low-temperature phase ^. We believe we have found important evidence of a doubly- 
degenerate ground-state. With the same density of states used to compute this evidence, 
we found this character changing as the temperature is increased towards the spin glass 
transition temperature. Just below this temperature, the behavior is consistent with the 
standard mean-field picture that has an infinitely degenerate ground state. 

The spin-glass state is characterized by randomly quenched exchange interactions with 
frustration. Often systems exhibiting this state are modeled by the ±J Ising model. In 
three dimensions this model is defined by the Hamiltonian 

n = -Y^ J^Jaiaj, (1) 

{id) 

where the Ising spins ai fluctuate thermodynamically, while each exchange interaction Jij is 
quenched to ±J randomly. The summation runs over the nearest-neighbor pair of sites on 
a cubic lattice. 

The frustration of local spin configurations by the exchange interactions generates many 
local minima in the free energy landscape, each minimum representing a seemingly random 
configuration of the spins. At low temperatures, the spins, however, may freeze into some 
configuration, and this freezing is the essence of the spin-glass "order." Because of the nature 
of this order, equilibration of spin glasses in simulations and experiments is often very hard. 

A. Nature of the spin-glass phase 

It is becoming increasingly accepted that the ±J model has a finite-temperature spin- 
glass phase transition in three dimensions, particularly after several recent Monte Carlo 
studies 0, 1^, 1^, ^. The nature of the low-temperature phase, however, is still controversial. 
Historically, the controversy has been mainly between advocates of the mean-field picture 
and the droplet picture. The mean-field advocates maintain the existence of an infinite 
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number of global minima of the free energy in the low-temperature phase. The condition is 
called replica symmetry breaking . The breaking of this symmetry is rigorously true for the 
mean-field (Sherrington-Kirkpatrick) model |Ty, |lT| of spin glasses. The question is whether 
it is also found in finite-dimensional non-mean-field models. The droplet advocates O, 



H| on the other hand, assert that the free-energy landscape has only two global minima 
that are connected through spin inversion symmetry. In the droplet picture, the nature of 
the three-dimensional ground state is seemingly less exotic than in the mean-field picture. 

The difference in the ground state degeneracy between the two pictures becomes quan- 
titative when we use the overlap order parameter of the spin-glass phase [|l|. To define 
this order parameter, we replicate the random exchange interactions in and change the 
Hamiltonian to 

E T.J^A"^<^f\ (2) 

a=l,2 (i j) 

where the superscript on the spin variables labels a replica. The overlap order parameter is 
then defined as the Hamming distance between spin configurations of the two replicas 0]: 



1 (1) (2) /„v 



i=l 

Here L is the linear size of the system and d is the dimensionality, which is three throughout 
the paper. 

The spins of one replica and those of the other are thermodynamically independent. Still 
they can be correlated because of the common random exchange interactions. In the high- 
temperature limit, the spin configurations are uncorrelated, and hence the overlap order 
parameter tends to zero as L~'^^'^. In the low-temperature limit, on the other hand, the spin 
configurations are frozen into energy minima, and hence the overlap order parameter can 
have a finite value. 

In the droplet picture, the overlap order parameter at low temperatures and in the ther- 
modynamic limit can take only two values equal and opposite to each other, as the frozen 
spin configuration of one replica is either macroscopically identical to the configuration of 
the other replica or to its inverted configuration. (We note, however, that the identity needs 
only to be macroscopic; the configurations can differ locally.) Hence the overlap order pa- 
rameter q takes one of the two values with an equal probability. In the mean-field picture, on 
the other hand, the overlap order parameter can take various values. Because of the many 
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free-energy minima, the frozen spin configurations of the two rephcas can be macroscopically 
different. 

More exphcitly, we define the order-parameter distribution as 



P{q) 



(4) 



^jD{E,q)e-^'^dE 
where the partition function is given by 

Z = Y, e-^^(M) = [dq [ dED{E, g)e-'^^. (5) 

{a} 

In these two equations D{E,q) is the normahzed density of states for the energy E and 
the overlap order parameter q. The energy in Eqs. (^) and and other energies hereafter 
(including those in the Monte Carlo simulations) are the energies of the replica Hamiltonians 
(0). The square brackets in Eq. denote the random average over various samples {Jij}- 
The distribution P{q) is normalized so that 

IPiQ)dq=^j:Piq) = l. (6) 
Physical quantities that are functions of q only are calculated as 

[(/(?))]av. = / mmdq, (7) 

where the angular brackets denote the thermodynamic average; e.g., the spin-glass suscep- 
tibility 

Xsg = L''f32 [{q2) - (g) 2]^^ (8) 

and the Binder parameter 

The droplet picture and the mean- field pictures have the functional forms of P{q) as shown 
in Fig. |l|. In the droplet picture, the order-parameter distribution in the thermodynamic 
limit has two delta-function peaks, indicating two possible thermodynamic states (Fig. |l]a). 
In the mean-field picture the order-parameter distribution is continuous between these peaks 
(Fig. |l|b). 

In practice, we can only simulate finite systems for which the order-parameter distribution 
looks like Fig. |l|c. We need to see whether the distribution in Fig. |l]c converges to Fig. ^ or 



U 
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FIG. 1: The functional form of the order-parameter distribution P{q) in the low-temperature phase: 
(a) According to the droplet picture; (b) According to the mean-field picture; (c) For finite-sized 
systems. 

Fig. |I|b as L — > oo. If linii^oo -P(O) — ^ 0, then significant doubt is cast upon the mean-field 
picture. 



While some simulations appear to support infinite degeneracy [0, [1^, ^7j, others do not 



?0|, and a recent analytic treatment argues against it |2^, We note 

that the droplet picture is a theory concerning the zero-temperature fixed point, while most 
studies suggesting the validity of the mean-field picture are based on numerical estimates of 
P(0) at T only as low as 0.7Tc, where the glass transition temperature is approximately 
1.0. If the doubly degenerate ground state predicted by droplet picture is correct, it should at 
least be seen at temperatures lower than 0.7Tc, the closer to T = the better. Unfortunately, 
the slow dynamics of spin glasses has hindered numerical simulations from exploring the 



vicinity of the zero-temperature fixed point. We comment that a recent numerical study ||25 
demonstrated that results of the Migdal-Kadanoff approximation appear to support the 
mean-field picture near and below the glass transition temperature, but eventually support 
the droplet picture as T ^ 0. Uncertainty, however, remains because the Migdal-Kadanoff 



approximation favors the droplet picture, even for the mean-field model [26, 27 



During the course of our work and just after, the possibility that the behavior of the spin 
glass does not unequivocally fit into one of these two pictures over the entire temperature 
range began to be appreciated. For example, Krzakala and Martin proposed the novel 
perspective that the entropy fluctuations lead to a trivial -P(O) at zero temperatures, even 
if there are zero-energy large-scale excitations (complex energy landscape). They further 
proposed that such a situation should arise in the the three-dimensional ± J Ising spin glass, 
and argued that if the energy landscape is complex with a finite number of ground state 
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families, then replica symmetry breaking reappears at finite temperatures. This perspective 
contested the standard picture. Palassini and Young ^ |2T[ studied this scenario and 
concluded that the size dependence of P{q) around g = is trivial and does not support the 
ultrametric picture. 

Katzgraber, Palassini, and Young pUl BD] addressed the issue of Gaussian versus bimodal 
distribution by simulating at finite temperatures the Gaussian exchange model and con- 
cluded that as in the bimodal exchange model -P(O) is trivial in the thermodynamic limit 
and suggested the existence of low finite energy excitations that cost finite energy and whose 
surface has fractal dimension less than the spatial dimension of the system. However, the 
sizes they studied were quite small and so they concluded there might be a crossover at 
larger sizes to a different behavior, such as a droplet or replica symmetry breaking picture. 

Two related papers by Hed et al. [0, Q suggest that the spin-glass phase possesses 



some characteristics of the mean-field description as a non-trivial P(0) and a hierarchical 
(but not ultrametric) structure of the pure states; nevertheless, they also claim that this 
phase is consistent with the Fisher-Huse scenario of the droplet picture. Correlated spin 
domains serve as the cores of zero enegy excitations. 



The recent work of Marinari et al. adopted the cluster analysis of Hed et al. ||3| 
and found strong continuity among physical features for T > and those found at T = 0, 
leading to a scenario with emerging mean-field like characteristics that are enhanced in the 
large volume limit for T > 0. These mean-field like features arise with entropic fluctuations. 



More recently Lamarcq et al. have studied the fractal dimension of the clusters that are 
the low-lying excitations of the model. 



There are still other papers, for example, ||36| , |37| , p8| , |39| , all illustrating a shift from 
simply the droplet versus the mean-field picture to something more subtle, with a consensus 
still very much evolving. One focus is on the nature and geometrical structure of ground- 
state excitations. 



B. Present Study 

In the present study, we significantly reduced the difficulty of the slow dynamics in Monte 
Carlo simulations by using a bivariate version of Berg and Neuhaus's multicanonical Monte 
Carlo method. Multicanonical simulations are performed independent of temperature or 



7 



of a range of nearby temperatures and estimate the density of states. From it we can in 
principle calculate expectation values at any temperatures. Our resulting estimates of -P(O) 
at low temperatures and as a function of lattice size suggest that at very low temperatures 
thermodynamic behavior could be consistent with a doubly degenerate ground state, while 
at higher temperatures, it is consistent with an infinite degenerate ground state. Hence, 
we found that an intrinsic temperature-independent quantity seemingly exhibited different 
looking equilibrium behavior at different temperatures. 

We present our Monte Carlo results in Sect. |T|. The order-parameter distribution function 
P(g) at the low temperature of T = 0.3 exhibits features indicative of a doubly degenerate 
ground state: P(0) decreases as the system size is increased. The low-temperature behavior 



of the Binder parameter also suggests double degeneracy. In Sect. |T| we discuss further 
implications of our results for the ground-state degeneracy. In Appendix, we describe our 
simulation method, namely a bivariate multicanonical Monte Carlo method. Monovariate 

we 



multicanonical methods [40, 41 1 have been applied to spin glasses before [42, 43, 44 



found, however, that by using a bivariate version we could reduce the correlation time 
of the simulation significantly. We show that the autocorrelation time is approximately 
proportional to the system size. 



II. NUMERICAL RESULTS 



We have carried out a bivariate multicanonical Monte Carlo simulation of cubic systems 
with edges L = 4 (1904 samples), L = 6 (2843 samples), L = 8 (1015 samples), and L = 10 
(1111 samples). The simulation method directly returns the density of states D{E,q), a 
temperature-independent quantity. We describe the details of the method and its use in 
Appendix; however, because the method is subtle, relatively new, and quite different from 
standard methods, we now summarize some of what is discussed in the Appendix. 

With the multicanonical method, sometimes called the entropic sampling method, we do 
not equilibrate the simulation at any value of T; that is, we never sample the steady-state 
distribution, generated by a Markov chain, that is supposed to represent the Boltzmann 
distribution. Instead, we sample from a steady-state distribution, generated by a Markov 
chain, that is adaptively constructed to be flat on the average. The flatness means we 
sample all accessible {E, q) values with an equal probability. In other words, we sample all 
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thermodynamics states equally. In fact, the sampling emphasizes regions where D{E,q) is 
small and as a consequence are generally difficult to access with many other methods. 

Sampling from this flat distribution also allows one to estimate D{E,q). We obtain 
D{E,q) for different system sizes. Once we obtain it for a given size, we can in principle 
obtain the the properties of the system for any temperature; that is, properties of the system 
at different temperature, no matter how different they may seem, all follow from the same 
D{E,q). 

The validity of our low T predictions depends on the "flatness" extending to the ground- 
state energy. As we discuss more fully in the Appendix flatness over the entire range of energy 
is not essential if not assumed in the evaluation of expectation values. The multicanonical 
method generally is a good ground-state sampler. This conclusion is supported in part by 
our to-be-reported rapid convergence of the entropy of the ground state for low T. As we 
will show, our estimates of the ground state residual entropy are consistent with the ones 



which Hartmann obtained by a ground-state counting algorithm |45]. This agreement is one 
reason why we are confident our simulations are covering low-energy states properly. At 
high temperatures, we will also see that the same D{E,q) estimates the Binder parameter 
to within two sigmas with the previous results 

A. Order-parameter distribution 

From the density of states D{E,q), we straightforwardly calculated the order-parameter 
distribution P{q), following Eq. (^). Figure ^ shows the temperature dependence of P{q) 
for L = 8. The function is close to a Gaussian distribution at high temperatures and has a 
double-peak structure at low temperatures. The results for T = 0.3, 0.4, and 0.5 are shown 
in Fig. 1^ for several values of L. We clearly see the decreasing tendency of -P(O) as L — > oo. 
The raw data, however, presumably underestimates the true values of -P(O): the P{q) data 
for L = 8 noticeably oscillates for g < 0.7 with the data point at g = happening to be in 
a valley of the oscillation. We presume that this oscillation is due to correlations between 
data at different values of q, for example, the data point P(0) being correlated with P(O.Ol); 



similar oscillations are seen in the data of others |2l|, qJ, ^ that was obtained by quite 
different methods. The present numerical method, as is summarized in Appendix, generates 
a random walk in the macroscopic phase space. The frequency of access by the random 
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€ 

FIG. 2: The temperature dependence of the order-parameter distribution P{q) for L = 8. The 
curves with the peak positions ordered from right to left correspond to T = 0.3,0.4,0.5, . . . ,2.0. 
The statistical errors are indicated only at each peak position. We plotted only for < g < 1, 
evoking the symmetry P(—q) = P{q). 

walker may be statistically less in some area of the phase space. The density of states will 
be underestimated in such area. Thus the correlations between data at different values of q 
can occur. 

To obtain proper estimates of -P(O), we have to smooth out the oscillations. We did this 
by choosing seven data points at the intervals of Aq ~ 0.1 over the range < g < 0.7 and 
then least-squares fitting them to the function InP(g) ~ cq + Cig2. The fitting parameter cq 
yields the estimate -P(O) = exp(co). The smoothed curves and the thus-estimated values of 
P(0) are shown in Fig. ^ for T = 0.3. We still see the decreasing tendency with increasing 
L. 

On the other hand, for T > 0.7, -P(O) does not show this tendency, even though its 
value is calculated with the use of the same D{E,q) (Fig. ||). In fact, -P(O) appears to 
converge to a finite value, a behavior which we now argue is spuriously consistent with the 
low-temperature behavior predicted by the mean-field picture. 

Independent of the degeneracy, we know from the scaling ansatz that at the critical point 
T = Tc ~ 1, -P(O) should increase as L oo. On the other hand, at very low temperatures, 
double degeneracy requires that -P(O) should decrease and infinite degeneracy requires that 
it should tend to a constant. What follows from our computed D{E, q) for different system 
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(e) (f) 
FIG. 3: The size dependence of the order-parameter distribution P{q). (a) A hnear plot for 
T = 0.5 and for L = 4, 6, 8, and 10. The peak position moves left as the system size is increased. 
Because the data points are very dense for L = 10, the error bar is shown only at the peak, where 
the statistical error is the largest, (b) A semi-logarithmic plot for T = 0.5. The error bars are 
shown only for a part of the data points, (c) A linear plot for T = 0.4 and for L = 4,6, and 8. The 
peak position moves left as the system size is increased, (d) A semi-logarithmic plot for T = 0.4. 
The error bars are shown only for a part of the data points, (e) A linear plot for T = 0.3 and for 
L = 4, 6, and 8. The peak position moves left as the system size is increased, (f) A semi-logarithmic 
plot for T = 0.3. The error bars are shown onl^-j^for ^ P^^^ of the data points. 
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FIG. 4: The fitting curves for P{q) and tlie estimates of P{0). Error bars are shown at only select 
points to illustrate the typical error in different ranges of q. 

sizes is shown in Fig. ^. Near T^, -P(O) does increase with an increasing L. At lower 
temperatures, however, it tends to a constant, a behavior supporting infinite degeneracy. At 
still lower temperatures, it tends to zero, a behavior consistent with double degeneracy. We 
note that a crossover scenario from critical, mean-field, to droplet behavior was predicted 



from the Migdal-Kadanoff approximation |25|. In particular we point out the similarity of 



our Fig. 6 to their Fig. 5. We thus suggest that most previous Monte Carlo studies, claiming 
to see behavior supporting the mean- field picture (infinite degeneracy), based on results for 
T only as low as 0.7, missed behavior consistent with the droplet picture (double degeneracy) 
which only appears at much lower temperatures. 



Palassini and Young have studied the scaling of the smoothed quantity x{l/2) 



1 /2 

J_i/2 P{(l)dq, and found a crossover scaling between T = behavior, where x(l/2) becomes 
trivial for L ^ oo and finite-temperature behavior, where the non-trivial part of P{q) has a 
much weaker dependence on L and is possibly size independent. The crossover is consistent 
with the qualitative features of our results. In fact, Palassini and Young's Fig. 4, showing 
this crossover, is strikingly similar to our Fig. |[ The remark by Palassini and Young that 
our P(0) drops dramatically at low T as L increase was based on the preliminary analysis 
of our data where P{0) was determined in the absence of the smoothing. Palassini and 
Young also point to possible different behavior between models with Gaussian distributed 
exchange interactions, where P{0) might be non-trivial as L — oo and bimodal ones, where 
P(0) becomes trivial. 



We also note the recent zero temperature work by Hed, Domany, and Hartmann [p4 . 
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FIG. 5: The size dependence of the order-parameter distribution P{q). (a) A hnear plot for 
T = 0.8 and L = 4, 6, 8, and 10. The peak position moves left as the system size is increased. 
Because the data points are very dense for L = 10, the error bar is shown only at the peak, where 
the statistical error is the largest, (b) A semi-logarithmic plot for T = 0.8. The error bars are 
shown only for a part of the data points, (c) A linear plot for T = 0.7 and for L = 4,6, 8, and 10. 
The peak position moves left as the system size is increased. For L = 10, the error bar is shown 
only at the peak, (d) A semi-logarithmic plot for T = 0.7. The error bars are shown only for a 
part of the data points, (e) A linear plot for T = 0.6 and for L = 4, 6, 8, and 10. The peak position 
moves left as the system size is increased. For ^3= 10, the error bar is shown only at the peak, (f) 
A semi-logarithmic plot for T = 0.6. The error bars are shown only for a part of the data points. 



FIG. 6: The size dependence of P{0) for T = 1.1, 1.0,0.9, . . . ,0.3. The data for T = 0.8 appear 
to be independent of L, but the data at lower temperatures reveal features of the droplet picture. 
The data in this figure were obtained after data processing ihustrated in Fig. ^ 



These investigators also found the need to smooth the values of P{q), and for better statistics 
they chose to study x* = J^ 'l P{q)dq. In the ground state they claim that x* scales to a 
small non-trivial value. To be more precise, they first separated P{q) into a part that comes 
from the big peaks close to g = ±1 that have an L-dependent tail at g = which scales to 
zero and into a part more proper to ground-state excitations whose scaling with L is the 
central issue. For this latter part they claim x* is non-trivial in the thermodynamic limit. 

B. Binder parameter 

The Binder parameter for spin glasses gsg, defined by Eq. (|), is essentially the kurtosis 
of the order-parameter distribution P{q). Because of the dimensionless combination of the 
second moment and the fourth moment, this parameter (except for effects due to correction 
to scaling) is expected to be independent of the system size at fixed points, i.e., T = 0, 
T = Tc, and T — oo . For conventional phase transitions such as ferromagnetic transitions, 
its temperature dependence for various system sizes has a crossing point at the critical tem- 
perature. At the high-temperature (T — > oo) fixed point, the order-parameter distribution 
P{q) should be Gaussian. Hence we should have 




(10) 
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FIG. 7: The temperature dependence of the Binder parameter g{T,L) for L = 4, 6, and 8. The 
statistical errors are comparable to the symbol size. 

or 

= 3 imi. . (11) 

Thus the high temperature fixed-point value of the Binder parameter is zero, i.e., 
(7sg(oo,L) = 0. In the high-temperature phase, the Binder parameter gsg{T,L) is renor- 
malized to be zero from above as L ^ oo. At the low-temperature (T = 0) fixed point, on 
the other hand, if the order parameter takes only two values ±go, as in the thermodynamic 
limit of usual ferromagnets, we have 









m 


av. 





and hence the fixed-point value of the Binder parameter is unity: Qsgi^, L) = 1. In the low- 
temperature phase, the Binder parameter gsg{T,L) is renormalized to be unity from below 
as L — > oo. At T = T^, the Binder parameter gsg{Tc,L) is expected to have a nontrivial 
universal value between zero and unity. Thus the crossing point of gsg{T, L) should give the 
critical temperature Tc. 

Our Monte Carlo simulation found the crossing point of the Binder parameter as shown 
in Fig. |4^. The critical point Tc should be in the region 0.8 >Tc> 1.1. For the moment, 
because of strong corrections to scaling, it is difficult for us to carry out sophisticated scaling 
analysis and obtain a more accurate estimate of Tc. Previous studies f^, H, ^ |^ on much 
larger systems claim Tc ~ 1.1. 
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FIG. 8: The temperature dependence of the Binder parameter g{T,L) for L = 4, 6, and 
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FIG. 9: The size dependence of the Binder parameter at T = 0.3. The Binder parameter is 
approaching unity rapidly. 

We now offer further evidence for a doubly degenerate ground state, as opposed to an 
infinitely degenerate one, by reporting our results for the low-temperature behavior of the 
Binder parameter in Fig. |^. These results strongly suggest that the Binder parameter tends 
to unity as T and L ^ oo. This behavior is consistent with Fig. |l]a and not with 
Fig. |l|b, as the Binder parameter is less than unity even at T = and L —>■ oo. Our Monte 
Carlo results thus clearly support double degeneracy (see Fig. 
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(a) (b) 
FIG. 10: The temperature dependence of the energy density and the free-energy density, (b) is an 
enhanced view of (a). The statistical errors are smaller or comparable to the symbol size. 

C. Residual entropy 

We calculated the entropy density from the difference between the energy and the free 
energy, i.e. 

(13) 

NT ' ^ ^ 

where F is the free energy calculated from the Monte Carlo output D{E, q) via 



F = -i log Z = -i log ( |: D{E, q)e-^^ ] . (14) 



Figure [1^ shows the energy and free-energy densities. The difference between them is shown 
in Fig. It is clear that at T = 0.1 the estimates of the entropy are virtually the resid- 
ual entropy at T = 0. The residual entropy is plotted in Fig. O. The convergence to the 



thermodynamic limit is quite rapid and the entropy seems to remain finite in the thermo- 
dynamic limit. Hartmann also observed a rapid convergence to a finite entropy from a 
ground-state search on systems up to L = 8. He obtained a similar estimate of the residual 
entropy s{T = 0) = 0.051(3)A;b. 

Although Hartmann used the existence of the residual entropy as evidence for the mean- 
field picture, its existence is on the contrary entirely consistent with the droplet picture. 
The degeneracy of the ground states predicted by these pictures is the degeneracy of ther- 
modynamic (macroscopic) states, while the residual entropy comes from the degeneracy of 
microscopic states. The distinction is important to note. 
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FIG. 11: The temperature dependence of the entropy density, (b) is an enhanced view of (a). The 
statistical errors are smaUer or comparable to the symbol size. 
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FIG. 12: The size dependence of the residual entropy. The statistical errors are smaller or compa- 
rable to the symbol size. 

Because the energy of a finite-sized ± J model is discrete, there is an inevitable degeneracy 
of the ground states. The issue is whether the degeneracy arises from these microscopically 
degenerate states or from many macroscopically different states. To make the distinction 
clearer, we will now consider a toy model in which we quench the exchange interactions 
into a periodic configuration with a unit cell of linear size /. This model has a ground state 
with a periodic spin configuration. We will assume, however, that every unit cell has one 
connected cluster of spins such that the spin inversion of the cluster does not change the 
ground-state energy (See Fig. 0). We refer to such spin clusters as "zero-energy droplets" 
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FIG. 13: Toy model with a finite residual entropy density and delta-function peaks in P{q). In the 
figure, L is the linear size of the whole system, / is the size of the unit cell, and v is the volume of 
a zero-energy droplet. 
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]. The number of such droplets in a system of linear size L is A^zed = (L/l) , where 
d is the dimensionality. The degeneracy of the ground-state energy is 2^^"='^. Therefore, the 



residual entropy density of the toy model takes a finite value 



Stoy(T = 0) =L-'^x A;Bln (2 



(15) 



This model, on the other hand, produces a P{q) consistent with the droplet picture: 
Consider two replicas of the toy model. Without loss of generality, we can fix the spin 
configuration of one replica and calculate contributions from different spin configurations of 
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the other rephca. The distribution of the order parameter is given by 

|g| = 1 with the probabihty 2"^^°^^ 

|g| = I — 2v/L'^ with the probabihty 2"^'"=^ x A^^ed, 

|g| = 1 - Av/L'^ with the probabihty 2"^-d x iiV^edlA^zod - 1) 



(16) 



2nv/L'^ with the probabihty 2 



A^zcd 

n 



where v is the volume of each zero-energy droplet. In the thermodynamic limit L ^ oo 
or equivalently A^^cd oo, the probability for n = A^zed/2 overwhelms and hence P{q) 
converges to a delta-function peak at 

q=l-^2vL' = l-^^. (17) 

and to a second delta-function peak at —q. 

III. GROUND-STATE DEGENERACY OF THE ±J MODEL 

The toy model just introduced is easily converted to a more realistic one for the 
ground-state degeneracy of the ±J model by allowing unequally-sized droplet volumes Vi 
{i = 1,2,3, ... , N^cd{L)). In defining the volume of these droplets, we always take the mini- 
mal volume; for example, if a zero-energy droplet contains another zero-energy droplet, we 



decompose them into two droplets as shown in Fig. |T^. In addition, the possible maximum 
volume of the zero-energy droplet is L'^/2: if there is a zero-energy droplet larger than the 
half of the system, we re-define the original ground-state spin configuration by fiipping the 
largest zero-energy droplet. 

The ground-state degeneracy is 2^^'='^'^'^^ and hence the residual entropy density is given 

by 

s{T = 0;L) = L-'^fcB In (2^-<^(^)) = L-'^iV,ed(^)A;B ln2. (18) 
We now define the distribution of the zero-energy droplets as rizediv) so that we have 

iV,ed(L) = L"^ n,,A{v)dv. (19) 
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• 



FIG. 14: If a zero-energy droplet contains another zero-energy droplet as shown on the left, we 
define the volumes of the two droplets as shown on the right. 



and rewrite Eq. ( ]I8| ) as 

s{T = 0; L) = kB\n2 X J n^ediv)dv. (20) 

Next, we consider the ground-state value of {\q\), and again we fix the spin configuration 
of one replica to the original ground-state spin configuration. Various contributions come 
from the spin configurations of the other replicas with some of the zero-energy droplets 
flipped. We extend Eq. ([T^) and obtain 

|g| = 1 with the probability 2^^='°'^, 

\q\ = 1 — 2vi/ L'^ with the probability 2~^^<=d for z = 1, 2, . . . , A^zed, 

|g| = 1 — 2(t>j + Vj)/ L'^' with the probability 2^^^°^ fQj- j = 1, 2, . . . , N^^d with i < j, 

|g| = 1 - {2/L'^) El^=i Vi^ with the probability 2-^-^ 

for ii, 22, • • • , «n = 1, 2, . . . , A^zed with ii < 12 < ■ ■ ■ < in, 



Summing all contributions, we find 

2 



(kl) = l 



2^zed Ijd 



i i<j 



1 — a 



zed 5 



where 



Q^zcd 



(21) 



(22) 



(23) 
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We can now compute (g2) in the same way and obtain 

2 1^ f 2 



(24) 



which reduces to 



or 



where 



Similarly we have 



(g2) = (l-a,ed)2 + /3,ed, (25) 



(g2)-(|g|)2 = /3,,d, (26) 



PzcdiL) = jj^Ylvi2 = — v2n^cdiv)dv. (27) 



where 



(g4)-3(g2)2 + 2(|g|)4=-27,ed, (28) 

lzed{L) = YmJ2^^^ = v4n^ed{v)dv. (29) 

If the order-parameter distribution P{q) has only two delta-function peaks at the L — > oo 
limit, then the left-hand sides of Eq. ( pB] ) and Eq. ( pH] ) vanish. This condition is satisfied if 
the density distribution of the zero-energy droplets n^cdiv) decays fast enough as v ^ oo. 
If n^cdiv) decays exponentially, the integrals in Eqs. ( jn\ ) and ( pPj ) give finite values; hence 
P^cdiL) = 0{L^'^) and 7zcd(-^) = 0{L~'^'^). If n^ediv) decays as with a; > 2, we have 
/3,od(^) = 0(L-'^'^''^(i'^-2)) 7,<,d(L) = 0(L-'^'^''^(3,a;-2)^)_ 

On the other hand, if the order-parameter distribution P{q) has a non-trivial part as 
L — > oo, azed(-Z^), /5zed(-Z^), and 7zcd(-Z^) all give finite values as L ^ oo. These conditions are 
satisfied if the distribution function decays as n^ediv) ~ f 

In our numerical results, both P^ediL) and 7zed(-^) appear to be decreasing as L is increased 
(Fig. |15D. The decrease is roughly 0{L^^), although the statistical errors are too large to 
make a more definite statement. This behavior is additional evidence for a doubly degenerate 
ground state. From our numerical results we suggest that the size distribution of the zero- 
energy droplets decays as rizediv) ~ with x > 2. 
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FIG. 15: The size dependence of (a) ^zediL) and (b) 7zed(-^)- Both quantities are vanishing as 
L — > CO. 



IV. SUMMARY 



The results of our multicanonical Monte Carlo calculation suggest a doubly degenerate 
ground state. Our main findings supporting double degeneracy are: (i) P{q) near g ~ 
decreasing at low temperatures as the system size is increased; (ii) the Binder parameter 
approaching unity at low temperatures; (iii) the effect of the ground-state degeneracy on 
moments of the overlap order parameter. In different temperature ranges, however, we 
further demonstrated, within finite size limitations, that the same density of states gave 
different scaling behaviors as a function of the system size. Only at very low temperatures 
is the scaling consistent with a doubly degenerate ground state. 

We believe our numerical results are consistent with recent analytic and numerical work 
and spin-glass themes pointing to a decoupling of the degeneracy of the ground state from 
the standard droplet and mean-field pictures and placing emphasis on the nature of the 
low-energy excitations and entropic fluctuations with them. Our numerical work provides 
information on these excitations only through the density of states. The P(0) predicted 
from this density of states is trivial at low temperatures. For finite systems, our density of 
states predicts crossover thermodynamic behavior as the temperature is lowered. 
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APPENDIX A: MULTICANONICAL MONTE CARLO METHOD 

We will describe the multicanonical Monte Carlo method from a viewpoint 0] slightly 
different from other authors [HOl O, H3, H3, M . We first compare it with random-sampling 



and canonical methods. 

Suppose that we perform a Monte Carlo simulation aimed at sampling from the weight 
W{E). After the simulation has equilibrated, that is, after the Markov chain reaches a 
steady-state, each microscopic state (spin configuration) is generated at a rate proportional 
to W{E). Also suppose that in "equilibrium" we construct a histogram h{E) of the energy 
E values associated with each configuration (state) generated. The rate at which a state 
with the energy E appears is proportional to D{E)W{E), where D{E) is the density of 
states. 

In a random-sampling simulation, one generates all microscopic states at the same rate. 
In other words, W{E) = const., and hence the resulting histogram is proportional to the 



density of states: h{E) oc D{E) (see Fig. 0). In many cases, the density of states becomes 
very small in low-energy regions. Hence simulations based on the random sampling generally 
have difficulty investigating low-temperature properties because the low-energy states, which 
dominate the thermodynamic average at low temperatures, are generated only infrequently, 
if at all. 

Importance sampling algorithms, like the Metropolis algorithm for canonical ensemble, 
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FIG. 16: Schematic of the histograms of energy values generated in the random-sampling, canonical, 
and multicanonical simulations. 

were developed to overcome this difficulty. Here states are generated with the importance 
weight W{E) oc e~^^ , and hence the energy histogram is h{E) oc e~^^D{E). Thus, in ideal 
cases, low-energy states appear at a greater frequency at low temperatures. In spin-glass 
simulations, however, the simulation tends to get stuck in local minima of the free-energy 
landscape. 

In the multicanonical simulations, the importance weight is set to W{E) oc 1/D{E), and 
hence the energy histogram should be fiat. Of course, one cannot set W{E) oc 1/D[E) 
because the density of states is a priori unknown. However, the multicanonical method is 
a procedure that makes the importance weight converge to the reciprocal of the density of 
states [^]. The aim is to ensure better statistics of low-energy states than in the random- 
sampling method and, at the same time, generate more high-energy states than in the 
canonical method so that the system can escape local free-energy minima. 

Berg and others have invented [^0|, ^ and applied P2| , this method to spin glasses, 



taking the weight of each spin configuration {a} either as 

W{{a})^l/D{E{{a})) (Al) 

or 

iy(M)cxl/D(g(M)), (A2) 

where D{E) is the density of states with respect to the energy and D{q) is the one with 
respect to the overlap order parameter. 
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In the ideal situation, the multicanonical simulation effectively generates a random walk 
in the macroscopic phase space. Because all the thermodynamic states should appear with 
the same probability, one can expect that the autocorrelation time of a thermodynamic 
variable {E in the above case) is proportional to the volume (the number of the sites). In 
the simulation of the ±J model, the energy change after each spin flip is of the order of 
J. On the other hand, the upper and lower bounds of the energy are of the order of NJ, 
where N = L'^ is the number of the spins. Hence the number of spin flips necessary for 
the random walk to cover the entire energy phase space is on the average of the order of 
A^^. Thus if random walk is in energy space, a one-dimensional space, the measured energy 
should decorrelate after A^^ Monte Carlo steps or N Monte Carlo sweeps. (A sweep is the 
attempt to flip each Ising spin once on average.) 

This ideal situation was not achieved in the previous applications of the multicanonical 



simulations to spin glasses [^2|, Berg et al. reported r oc A^ for the importance 



weight (^) 1 42 1 and r oc A^^-^^ for the importance weight (^) [44|. (Here r is measured 



in the unit of Monte Carlo sweep.) This suggests that the slow dynamics in the low-energy 
region was not removed completely in these mono-variate multicanonical simulations. (Note, 
however, that Berg et al. did not use the conventional definition of the autocorrelation time 
in measurement of r. This might be a contribution to the above power law behavior.) 
In the present study, we carried out the simulation using the bivariate importance weight: 

Wo^l/D{E,q). (A3) 

As shown in Fig. |1^ (a), we almost achieved the ideal situation t oc N. The random walk 
in the two-dimensional phase space is illustrated in Fig. |I7| (b) . The autocorrelation plotted 
in Fig. (a) was computed as follows: We first measured the autocorrelation as the simple 
Monte Carlo average without any weighting: 

CE{t) ^ -^EmEU + t) (A4) 

C,{t) ^ j^E^iMj+t), (A5) 

where A^mcs is the number of Monte Carlo steps for measurement. Then the autocorrelation 
time was computed as the integrated time: 

^E^^^ECEit), (A6) 

^E[yJ> t=0 
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FIG. 17: (a) System-size dependence of the autocorrelation time r with respect to E and q. The 
dotted Une represents the power law r oc A^^'^^, which was reported in (b) An example of the 
random walk in the phase space {E, q). The system size simulated here is L = 4. The time series 
of 2500 Monte Carlo sweeps after equilibration is plotted. 



r, 



' CM to 



(A7) 



where Ie (ig) is defined as the time when CE(t) (Cg(t)) first becomes negative. 

There is a drawback with the bivariate simulation: large electronic memory is needed to 
store the data from which the bivariate histogram is constructed. For the Ising model we 
used the number of all possible discrete states to equal the number of bins in the energy 
phase space. This number was of the order of L'^ and was similarly sized in the q phase space. 
Hence an array of the size L^'^ was necessary to store each density of states D{E,q), the 
histogram h{E,q), and related work space. The combined storage requirements amounted 
to nearly 30MB per sample for L = 10 and 400MB per sample for L = 16. Based on 
direct comparisons with the original mono-variate multicanonical method, we have, however, 
concluded that for spin-glass simulations the use of both E and q is essential for improving 
the slow dynamics 13] . We comment that bivariate multicanonical Monte methods have been 



previously used in different contexts, e.g., protein folding simulations |£0|, ^ 



The actual algorithm of the multicanonical method has the following structure: 

(A) Make a guess at the density of states, Dq{E, q). Start the following loop with z = 0. 

(B) The outer loop (multicanonical- iteration loop): 
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Set the importance weight to Wi = 1/Di. Calculate the energy and the order 
parameter of the initial spin configuration, Eq and go- Start the following loop 
with j = 0. 

The inner loop (spin-update loop): 

(i) Choose a spin to be updated. Calculate the energy and the order parameter, 
E'j and q'j, assuming that the spin is actually fiipped. 

(ii) Calculate the spin-fiip probability as 

V Wi{Ej,q,)J 

(iii) Draw a random number i? G [0, 1]. If -R < PfHp, fiip the spin and make the 
substitution Ej^i = Ej and g^+i = q'j. Otherwise, keep the current spin 
configuration and make the substitution Ej+i = Ej and qj+i = qj. 

(iv) Increment the histogram bin hi{Ej^i,qj^i) by one. This step should be 
skipped during the burn-in stage (specifically given below), that is, until the 
simulation reaches "equilibrium" for the importance weight W^. 

(v) Increment the number of steps j by one and go to the step (i) until j exceeds 
a pre-determined number. 

Measure physical quantities Q{E,q) via 

{Q)^ = ^ E QiE, g)e-^^/i.(i?, q)D,{E, g), (A9) 

E,q 

where the partition-function estimate is given by 

Z^ = Y. e-^'^HE, q)D,{E, q). (AlO) 

E,q 

In practice, this step was skipped until the histogram hi became reasonably flat 
(see below for more specifics). 

Update the guess of the density of states as follows: [E2[ 



D<E,q) = A(i^,g) ME,q) + l) ^ 
* ' Normalizaion 

The normalization constant is chosen so that the integration of the density of 

states becomes unity. (An extra count is added to the histogram in (|A11 ) in order 



to avoid zero division in setting the next importance weight VFj+i = l/-Di+i [0.) 
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(e) Increment the number of multicanonical iterations i by one and go to the step (a), 
until the histogram meets certain conditions of being flat. 

The method is an iterative procedure requiring an initial guess Do{E, q) of D{E, q). We 
found our average results to be quite independent of this guess, but used the following 
particularly convenient and efficient choice: Our choice was a flat one for L = 4, a case 
where the convergence was very rapid. The results for this system size were then scaled and 
used for the starting points for larger systems. 

For a given iteration, step (iv) should be skipped until the simulation reaches "equilib- 
rium" for the importance weight Wi. For L = 4, 6, 8, and 10, we used 10 000, 40 000, 50 000, 
and 100 000 Monte Carlo sweeps for this burn-in stage. The number of the burn-in sweeps 
were approximately a hundred times longer than the autocorrelation time given in Fig. 17a. 
The autocorrelation times were determined before we began our production running. After 
the burn-in we used 1 000 000, 4 000 000, 5 000 000, and 10 000 000 sweeps for the different 
systems sizes during each iteration step; that is, the number of sweeps in each iteration was 
10 000 times the measured autocorrelation time and one hundred times the burn-in time. 

Step (c) was skipped until the histogram hi became reasonably flat. The number of 
iterations required for this varied as function of L and from one random bond configuration 
to another for a given L. For L = 4, the range was 2 to 15; for L = 6, 2 to 23; for 
L = 8, 2 to 276; and for L = 10, 4 to 227. True flatness is never achieved nor is it 
necessary. What is important is the histogram be reasonably flat, extends to the band 
edges, and lacks spikes at the band edges. In defining a "reasonable" fiatness criteria, which 
is the criterion for the iteration to have converged, we were motivated by the following 
observations: While converging, the histogram typically is first concentrated in the mid- 
part of the density of states and expands outward towards the upper and lower bounds 
of the energy (band edges) and the order parameter until it becomes flat. At the leading 
edges, spiky structures sometimes appear as the histogram is crossing over from occupied 
to unoccupied states. When the edges are reached, the spikes disappear. We defined the 
histogram to be flat, whenever the total count at the minimum energy and that at = 
were within one sigma of the average histogram count. In contrast to the usual Ising model, 
the band width for a bond configuration is not known a priori, but it is easy to estimate 
from short runs and then set it slightly larger than expected to handle extremal states. In 
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general, we found the multicanonical method to be an efficient estimator of the lowest energy 
of a given bond configuration. 

After convergence, we repeated the outer loop five times after so that we could detect 
erroneous estimates of the density of states; that is, using the just determine D{E,q) as 
the starting point, we repeated multicanonical procedure until five consecutive iterations 
maintained convergence. Testing determined that five was an adequate number to avoid 
false convergence. We also divided the inner loop into ten blocks, from which we calculated 
the statistical errors. Because the length of each block is about 1000 times longer than the 
auto-correlation time, the absence of correlations between block averages is expected and 
standard estimates for statistical error were used. 

Because it is an importance estimator, ([A9|) could in principle be applicable before the 
histogram becomes flat; that is, achieving flatness is not an essential requirement for the 
accuracy of the method if importance weighted estimators are used. Even with convergence, 
we used such estimators for the sake of caution. We only collected data for converged results. 
With sufficient statistics, the histogram hi would be proportional to Wi{E,q)D{E,q) = 
D{E,q)/Di{E,q), where D{E,q) is the true density of states. Thus the estimator ( |A9|) 
would give 

{Q)^ ^ 1 5: QiE, q)e~^''D{E, q), (A12) 

E,q 

with the partition function 

Z = ^e-^^D(E,g). (A13) 

E,q 

We comment again that we could have calculated all physical quantities after the simu- 
lations were completed, provided we had stored the density of states for each random bond 
configuration on disk. This procedure would requires a huge amount of disk space because 
of the need to simulate a large number of random samples. 



We also note that Marinari et al. ||53| have recently tried to implement our bivariate ver- 
sion of Berg and Neuhaus's multicanonical method and reported experiences very different 
from ours. They reported difficulties in ergodicly sampling the importance function, finding 
it only remained "fiat" only by using an order of magnitude larger number of Monte Carlo 
sweeps than we used. This would correspond needing burn-in times and block sizes 1000 
times our estimated auto-correlation time. More importantly the very small values of hi er- 
ratically riddling their histograms would translate into very small values of hi{E, q)Di{E, q), 



30 



which should have caused our measured values for a given bond configuration to vary widely. 
We did not observe this, and our estimates of autocorrelation time were consistent for dif- 
ferent system sizes. As mentioned above, most of our data analysis was done "on-the-fly" 
to avoid storage problems so we cannot re-examine the actual data for our reported results. 
However, when a referee pointed out to us the Marinari et al. difficulties, we immediately 
executed a number of simulations for the L = 8 system size causing Marinari et al. prob- 
lems and looked at the same quantity that they reported in their Fig. 3. We never saw the 
behavior that they claimed occurs. The origin of the discrepancy between their results and 
ours is unknown to us. 
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